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Abstract 

This paper presents a novel mesh-update technique for unsteady free-surface Newto- 
nian flows using spectral element method and relying on the arbitrary Lagrangian- 
Eulerian kinematic description for moving the grid. Selected results showing com- 
patibility of this mesh-update technique with spectral element method are given. 
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1 Introduction 



Incompressible free-surface flows are encountered in a wide range of engineer- 
ing and environmental flows. In the nineties the more specific case of turbulent 
free-surface flows started to be investigated with numerical computation based 
on high-order methods [8,9]. In our work, we aim at computing large-eddy 
simulation (LES) of unsteady, incompressible and Newtonian turbulent free- 
surface flows by using the spectral element method (SEM) [13,14]. The choice 
of interface-tracking technique was made to ensure an accurate description of 
the free surface. 



This paper highlights the computational techniques we are developing for sim- 
ulating incompressible free-surface flows using the SEM. These techniques in- 
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elude the arbitrary Lagrangian-Eulerian (ALE) formulation [2,7,15], mesh 
update and re-meshing methods [6,10]. 

This paper is organized as follows. The governing equations in the ALE frame- 
work for general free-surface flows are introduced in Section [2j Then, we 
present the discretization methods and numerical technique in Section [3j Sec- 
tions H] and [5] are dedicated to the moving- grid technique and the mesh-transfer 
operation, respectively. 



2 Governing Equations 



A moving boundary-fitted grid technique has been chosen to simulate the free 
surface in our computations. This choice of a surface-tracking technique is 
primarily based on accuracy requirements. With this group of techniques, the 
grid is configured to conform to the shape of the interface, and thus adapts 
continually (at each time step) to it and therefore provides an accurate de- 
scription of the free surface to express the related kinematic and dynamic 
boundary conditions. 

The free-surface incompressible Newtonian flows that we have considered are 
governed by the Navier-Stokes equations comprising the momentum equa- 
tion and the divergence-free condition. In the arbitrary Lag- rangian-Eulerian 
(ALE), a mixed kinematic description is employed: Lagrangian description of 
the free surface c?f2p(£), Eulerian description of the fixed domain boundaries 
dVlo and mixed description of the internal fluid domain Q(t), subset of M. d 
with d — 2, 3 the space dimension, t referring to the time as the fluid domain 
is changing when its boundaries are moving. Let us denote by Qq a reference 
configuration (for instance the domain configuration at initial time t = t ). 
The system evolution is studied in the time interval I — [to,T]. The position of 
a point in the current fluid domain Q(t) is denoted by x (Eulerian coordinate) 
and in the reference frame f2 by Y (ALE coordinate). Let A be a family of 
mappings, which at each t E I associates a point Y e Q to a point x G Qf- 

A t ifloCR^-tfliC R d , x(Y,i)=A(Y). (1) 

At is assumed to be continuous and invertible on Q and differentiable almost 
everywhere in /. The inverse of the mapping At is also continuous on fl . With 
these notations the set of equation reads: 



<9v 
~dt 



+ (v - w) • V x v = -V x p + 2z/V x • D x (v) + f in Q{t), (2) 

Y 

V x -v = in Sl(t), (3) 
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with v(x, t) the velocity field, p(x, t) the pressure field (normalized by the 
constant density p), D x (v) = |(V x v + V x v T ) the rate-of-deformation tensor, 
v the kinematic viscosity of the fluid and f the body force. The ALE mesh 
velocity w(x, t) appearing in ([2]) is defined as 



w(x,t) = - 



dA, 



dt 



(4) 



Surface tension effects are assumed to be negligible as we deal with turbulent 
flows. The associated boundary conditions are: 

— the kinematic boundary condition on dVLp{t): 

v • n = w • n, (5) 

n being the local outward unit normal to the free surface; 

— the dynamic boundary condition on dflp(t): 

-pn + 2z/D x (v) • n = 0, (6) 

assuming an inviscid air and zero ambient pressure; 

— homogeneous Dirichlet boundary condition on dVtp,: 

v = w = 0. (7) 

In addition to the set of governing equations d^D - ©, the closure of this free- 
surface problem based on a moving-grid formulation requires one more equa- 
tion governing the evolution of the mesh velocity w in the internal fluid do- 
main Q(t). The boundary values of w being prescribed by the equations <^ 
and ([7]) on the boundary d$}p(t)Udflp> of the fluid domain. This last governing 
equation for w will be presented in detail in Section HI 

As our focus is on transient problems, proper initial conditions at time t = to 
for the fluid velocity v and for the mesh velocity w have to be provided. The 
initial fluid velocity must satisfy the divergence-free condition and the values 
of the initial mesh velocity have to be given together with the initial shape of 
the free surface. 

Based on the strong formulation of this free-surface problem given above, one 

can derive the more appropriate weak transient ALE formulation: 

Find (v(t),p(t)) G Hq D (tt(t)) d x L 2 (f2(t)) such that for almost every t > t 



— [ (u o at 1 ) ■ vdn + / (u o at 1 ) ■ v x • [w - vw] an = 

at Ju(t) ' Jn(t) 

[ ( P V X ■ (u o AT 1 ) - 2//D x ((u o AT 1 )) : V x v) dfi (8) 

Jn(t) 

+ [ f -{uo AT 1 ) aft WueH^ D (Q ) d , 
Jn(t) 
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and 



-/ gV x -vdfi = \/q eL 2 (n(t)). (9) 
Jn(t) 

with the functional space HQ D (Q(t)) defined by 

It is worth noting that the weak formulation (jSJ)-© is only valid in our par- 
ticular case where homogeneous natural and essential boundary conditions, 
respectively (jSJ) and (j7]) are applied to the system. 



3 Numerical technique and discretization 



A classical Galerkin approximation is applied to the set of governing equations 
in its weak transient ALE form (JHJ) on the flow domain Q(t), in order to 
determine the pressure and the fluid velocity, keeping in mind that the mesh 
velocity is obtained by the moving-grid technique developed in the next sec- 
tion. The Galerkin approximation is then discretized by using the spectral 
element method with the classical staggered P^r — Ptv-2 approach to avoid 
the development of spurious pressure modes. Discontinuous and continuous 
approximations are respectively taken for the pressure and fluid velocity. The 
mesh velocity is discretized using the same polynomial space as the fluid veloc- 
ity, namely Pat, based on a Gauss-Lobatto-Legendre (GLL) grid of order N. 
For the discontinuous approximation of the pressure, a Gauss-Legendre (GL) 
grid of order iV — 2 is used. Consequently the ALE Navier-Stokes semi-discrete 
equations can be derived from (jHJ)-©: 

— (Mv) + C(v,w)v = -Kv + D T p + F, (10) 
at - 

-Dv = 0, (11) 

M denoting the mass matrix, K the direct stiffness matrix, D T the discrete 
gradient operator, D the discrete divergence operator, C(v, w) the discrete 
convective operator depending both on the fluid and mesh velocities and F 
the discrete body force. The update of the position x of the mesh points is 
performed by integrating the following discrete equation: 

dx , . 

(12) 

The set of semi-discrete equations (|T0|) - (|T2|) is discretized in time using a 
decoupled approach: the linear Stokes computation (linear viscous diffusive 
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term) is integrated based on an implicit backward differentiation formula of 
order 2, the nonlinear convective term is integrated based on a simple method 
used by Karniadakis et al. [11], consisting in an explicit extrapolation of order 
2. Finally the update of the position of mesh points is based on an explicit 
and conditionally stable Adams-Bashforth of order 3 (AB3). 

Lastly the treatment of the pressure relies on a generalized block LU decom- 
position, using a standard fractional-step method with pressure correction. 



4 Moving-grid technique 

As already mentioned in the previous sections, our free-surface flow compu- 
tations are of interface-tracking type and rely on a moving-grid technique, 
allowing large amplitude motions of the free surface, generating a grid con- 
forming to the shape of the free surface for an accurate and easy application 
of the boundary conditions on <9fii?(t). Moreover a description as accurate as 
possible of the turbulent free-surface boundary layer is essential to our work. 
These points justify by themselves the choice of a moving-grid technique that 
increases the difficulty of the marginally intractable problem of turbulent vis- 
cous flow computations. 

The computation of the mesh velocity w in the internal fluid domain Q(t) is 
the corner-stone of the moving-grid technique developed in the framework of 
the ALE formulation. The values of the mesh velocity being prescribed on the 
boundary dfl(t) = dQp(t) U OVLd as expressed by equations ^ and ([7]), the 
evaluation of w in Q(t) can be obtained as the solution of an elliptic equation: 

£ x w = in tt(t). (13) 

This elliptic equation constitutes a classical choice for calculating the mesh 
velocity [8]. In the present case it is desirable to impose an additional con- 
straint to the mesh velocity problem, in order to ensure the incompressibility 
of the mesh by imposing a divergence-free condition to w: 

V x -w = in Q(t). (14) 

Our choice for the elliptic operator £ is based on the assumption that the 
motion of the mesh nodes is equivalent to a steady Stokes flow, corresponding 
physically to an incompressible and elastic motion of the mesh. The boundary- 
value steady Stokes problem for the mesh velocity can be formulated as follows: 

w ■ n = v ■ n on dVLp(t), (15) 

w-r = ondtt F (t), (16) 

w = on dn D , (17) 
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where r is the local unit vector directly orthogonal to n, and 



V x • cr = 
V x ■ w = 



in Q(t), 
in Q(t), 



(18) 
(19) 



denoting by cr the Cauchy stress tensor of the mesh defined by 



<T = — 



■pi + v{ V x w + V x w T ) 



(20) 



with p and v being respectively the fictitious mesh pressure and the fictitious 
kinematic viscosity of the mesh, characterizing the elasticity of the mesh in 
its motion. 

The choice of this boundary-value problem for the mesh velocity has several 
justifications. Constraining the elliptic equation by a divergence- free condi- 
tion for w allows to ensure the conservation of the volume of the spectral 
elements, condition that is helpful in practice to have rapidly convergent com- 
putations [1]. In general the global volume of the computational domain may 
not be conserved, e.g. with an inflow-outflow imbalance, which requires (fITl) to 
be relaxed. In addition, the mesh velocity w appears in the convective part of 
equations §2$), (jHD and (ITU]) , together with the divergence-free fluid velocity v. 
Moreover it is worth remembering that the divergence-free condition imposed 
to w leads to a conservation of the metrics (the Jacobian being constant in 
time) when moving the mesh. Finally the unavoidable issue of fulfilling the ge- 
ometric conservation law (GCL) in the ALE framework [3-5] is automatically 
solved when considering a divergence-free mesh velocity as a consequence of 
the work of Formaggia and Nobile in [5] . 

From a numerical point of view, the problem corresponding to the set of 
equations (IT5"|) - (IT9"|) is discretized using the SEM, with a staggered grid Pa? — 
Pat_2 for the couple mesh (w, p). An Uzawa decoupling technique is employed 
for the treatment of the fictitious pressure. 

Based on the technique described earlier, we have developed the following 
moving-grid algorithm: 

(1) Input data: mesh M. n aXt = t n , with nodal coordinates x n , fluid velocity 
v n on dO,p, mesh velocity w n in Q n U dVt n ] 

(2) Step 1: steady Stokes computation of w n+1 by Eqs. f|T5|) - f|T9|) : 

(3) Step 2: update of the nodal coordinates Eq. ffl2l ; spectral element vertices 
are moved according to the AB3 scheme: 



(4) Creation of the new mesh M. n+1 with the new Gauss-Lobatto and Gauss- 
Lobatto-Legendre grids for each new spectral element; 




(21) 
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(5) Output data: mesh M. n+1 at time-step t = t n+ i, with nodes coordinates 
x™ +1 , mesh velocity w" +1 in fi n+1 U dtt n+1 . 

Two performance tests have been carried out on a study case where one edge 
of a squared mesh is deformed by a sine profile. Both of these tests aimed 
at verifying the spectral element volume conservation that is theoretically 
imposed by the divergence-free condition on w. The first test is dedicated to 
the verification of the global volume conservation, by computing the relative 
change of the volume of the computational domain when moving the grid from 
the initial square to the deformed one. For several number of spectral elements 
and for a polynomial interpolation order ranging from 1 up to 12, the relative 
change of the volume of the computational domain is found to be smaller 
then the machine precision. The second test is also devoted to the volume 
conservation but now from a local perspective and by numerically computing 
the L 2 (Q)- and L 2 (w)-norm of the divergence of the mesh velocity w for a 
polynomial interpolation order N ranging from 5 up to 12, where u is interior 
of the computational domain made of the spectral of elements of Q not sharing 
an edge with dfl. Results are presented on Fig. [1] and it is found that these 
norms are exponentially decreasing with N as expected when using a spectral 
element method [1]. Moreover we can note that the L 2 (u;)-norm of V • w has 
a faster rate of convergence than the L 2 (fi) -norm. This is justified by the fact 
that the divergence-free constraint cannot easily be enforced at the grid points 
located in the vicinity of the boundaries of the computational domain Q. 



le-02 




le-06 r 

le-07 r ' •+.. '- 

le-08 - ' "■+..... 

'"■■+-... 

le-09 r " + 

le-10 t 1 1 1 1 1 1 1 

5 6 7 8 9 10 11 12 13 

N 

Fig. 1. L 2 -norms of the divergence of the mesh velocity w versus polynomial inter- 
polation order N (Log scale) 
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5 Mesh-transfer operation 



In the previous section was presented the moving-grid technique used in our 
work to move the grid points at each time-step, generating a new mesh. De- 
pending on the amplitude of the mesh deformation at each time-step, this 
technique can be applied during an important number of iterations. Never- 
theless the mesh obtained by moving the grid nodes can be too convoluted 
therefore affecting the accuracy and the convergence of the simulation. Conse- 
quently a re-meshing operation is to be called by a specific control parameter 
(e.g. a discrete Jacobian positiveness criterion) to provide a new mesh topol- 
ogy. Before starting the ALE Navier-Stokes computation at the next time-step 
on this newly created mesh, it is mandatory to transfer some information from 
the previous mesh to the new one. The main requirement imposed to this so- 
called mesh-transfer operation is to conserve the spectral accuracy of the SEM. 
The information to be transferred comprises six fields: the fluid velocities v n , 
v n_1 and the mesh velocities w n , w n_1 , w n ~ 2 (time-integration schemes are 
of order 2 for v and 3 for w) and also the pressure at the current time-step 
(use of a pressure correction technique). As written in Section [31 the velocities 
are expanded over a GLL grid and the pressure over a GL one. Therefore our 
mesh-transfer technique must be capable of transferring fields defined over GL 
and GLL grids. 

Our mesh-transfer algorithm for GL grids being based on the one for GLL 
grids, we will start presenting in detail the latter. Let us consider two meshes 
At 1 and Ai 2 corresponding to different mesh topology of the same computa- 
tional domain and the mesh-transfer operation from Ai 1 to At 2 . In the sequel 
we will assume that we have the following decompositions in terms of spectral 
elements: 

Ei 

QiUdQi^ U°*' e for z = 1,2. (22) 

e=l 

As the computational domain remains unchanged, for each spectral element 
f2 2,e of Ai 2 we have: 

^cffiiU^) Ve = l,...,£ 2 . (23) 

Due to Eq. (|2"3"|) our mesh-transfer technique only requires an interpolation 
procedure. Let us note the physical location of the set of GLL grid points of 
a spectral element Q 2,e2 (e 2 = 1, ■ ■ ■ ■ , E 2 ) by {x 2 j 62 } with (i = 1, ■ ■ • , N X)2 + 
l;j = 1, ■ • • , N y> 2 + 1), N x> 2 (resp. N V)2 ) being the order of the polynomial 
interpolation in the x-direction (resp. ^/-direction) for the mesh Ai 2 (with the 
same notations, N Xt2 and N y> 2 can be different from N x> i and N y x respectively). 
The proposed algorithm can be summarized in three steps: 

(1) Find the spectral element Q 1,ei of At 1 containing x 2 j 62 ; 
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1,61 of x^' 62 within the parent element Cl l,ei of 



(3) Compute the value of the field at the point x?j 62 given r 1,ei , the GLL 
Lagrangian interpolation basis and the values of the field at the GLL 
grid points of Q} ,ex . 

The first step causes no difficulty in its implementation. The second step uses a 
transfinite interpolation procedure in each spectral element, in order to invert 
the iso-parametric mapping 



r l,ei = ( r l,e-l ; s l, ei j = $-l( x 2,e 2 ) w j th pl.ei £ ^l.ej = ^ ^ ^ 



In practice, the inversion is carried out differently depending on the topology 
of the spectral element. With quadrangular spectral elements, our algorithm 
performs a direct analytical inversion of the affine mapping <fr which is compu- 
tationally inexpensive. With deformed spectral elements [1], the inversion of 
<]? relies on the so-called 'inverse iso-parametric mapping technique' from Lee 
and Bathe [12] which is based on a Newton-Raphson type iterative procedure. 

Finally in the last step, efficient routines compute the following spectral in- 
terpolation: 



with {7Tj(0}j=o an d P — x,y-, the one-dimensional GLL Lagrangian inter- 
polation basis of degree N Pt \. As said earlier the mesh-transfer technique for 
GL grids relies on the one for the GLL grids. In our simulations, the only 
GL-interpolated field that has to be mesh-transferred is the pressure field. 
Therefore, by interpolating the pressure on the GLL grid, then by applying 
the GLL mesh-transfer operation introduced earlier and finally by interpolat- 
ing back on the GL grid, we manage to perform the requested operation. It 
is important to minimize the occurrence of a re-meshing as our mesh-transfer 
technique is computationally expensive even for quadrangular elements (affine 
iso-parametric mapping). A more detailed assessment of the performance of 
this technique is provided at the end of this section. 

This mesh-transfer operation has been extensively tested in order to ensure 
its compatibility with the SEM, regarding its exponential rate of convergence. 
Tests involving the following two key parameters have been carried out: the 
polynomial interpolation order N and the amplitude of the change in topology 
of the grid when re-meshing. 

The set-up is presented in Fig. [2] and is made of a mesh comprising four spectral 
elements. The change in topology of the mesh is prescribed by moving only 



u ( r l.ei >s l,ei) =J2J2 U^r 1 ^)^ 1 ' 61 ) 



(25) 



k=0 1=0 
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the vertex u (see Fig. [2]) common to all four spectral elements and afterwards 
the mesh-transfer operation is performed. 




Fig. 2. Sketch of the computational domain O, the two meshes M 1 and M 2 and their 
spectral element decompositions before and after a prescribed re-meshing operation 
obtained by moving the central vertex oj 



To evaluate the dependence of our technique with the interpolation order N, 
the central vertex is moved to produce a topological change in the mesh by 
a factor of approximately 10 %. An analytical field / is interpolated on the 
initial mesh and mesh-transferred onto the distorted mesh, leading to the 
interpolated field /. The interpolation error is defined by e = \\f — f\\L 2 (n) 
and computed values are presented in Table [lj showing a conservation of the 
exponential rate of convergence. 
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N 


e = 


11/ - /IU 2 (n) 
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e = 


11/ - f\\iP(a) 
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Table 1 

Evolution of the error e with the spectral interpolation order N 



To characterize the effect of the distortion of the mesh on our mesh-transfer 
operation, all possible positions of the moving vertex within the computational 
domain Q were envisaged. In particular, we present here the case where uj is 
moved along the diagonal AC of the computational domain Q as shown in 
Fig. [2J Its motion is characterized by the set of coordinates (a, (3) of u in the 
parent domain Cl = [—1, l] 2 . The interpolation error e was again computed for 
three values of iV and results appearing in Table [2j show that our technique 
is totally independent on the amplitude of topological change of the mesh due 
to the re-meshing operation. 
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Table 2 

Evolution of the error e when uj moves along the diagonal AC, with 7 = 1 and for 
three different values of N 

Lastly, the computational expense of the mesh-transfer has been evaluated for 
a polynomial degree N = 10 in both directions (19 2 grid points for this 2D 
grid), and as previously for a topological change in the mesh by a factor of ap- 
proximately 10 %, corresponding to a 'small' 2D case. The results confirm the 
afore- mentioned cost: a complete mesh-transfer corresponds to approximately 
100 Navier-Stokes solves depending on the value of the time-step. 
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6 Conclusion and future studies 

A novel isochoric moving-grid technique and mesh-transfer technique for spec- 
tral element grids have been presented. Both of these techniques are the corner- 
stones of our computations of turbulent free-surface flows using spectral ele- 
ment method. Part of the work was to ensure that these two techniques have 
no effect on the exponential rate of convergence, the main reason of our choice 
of the spectral element method. We have obtained positive results all along 
the extensive series of tests carried out to verify the behaviour of this rate of 
convergence. The development of an automatized re-meshing scheme coupled 
to a re-meshing control parameter are still under investigations. 

Our next goal is to simulate three-dimensional turbulent free-surface flows 
using the techniques presented in this paper with the difficult task of gaining 
a better insight into the physics involved in the thin turbulent boundary layer 
near the free surface. 
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